Optimal waveform for the entrainment of a weakly forced oscillator 
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A theory for obtaining waveform for the effective entrainment of a weakly forced oscillator is 
presented. Phase model analysis is combined with calculus of variation to derive a waveform with 
which entrainment of an oscillator is achieved with minimum power forcing signal. Optimal wave- 
forms are calculated from the phase response curve and a solution to a balancing condition. The 
theory is tested in chemical entrainment experiments in which oscillations close to and further away 
from a Hopf bifurcation exhibited sinusoidal and higher harmonic nontrivial optimal waveforms, 
respectively. 
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Entrainment of oscillators to an external signal in non- 
linear dissipative systems is a fundamental concept of im- 
portance in a large variety of applications [1]; two promi- 
nent examples include the time-scale adjustment of cir- 
cadian system to light [2] and the cardiac system to a 
pacemaker [3| . The entrainment process is rigorously de- 
scribed theoretically by phase/amplitude equations and 
circle maps for weakly and strongly perturbed nonlinear 
systems, respectively [1]. The general result of the the- 
oretical analysis is that nonlinear oscillators can adjust 
their frequencies to that of the external source above a 
critical forcing amplitude. In the forcing amplitude vs. 
forcing frequency diagram there are long vertical entrain- 
ment regions called Arnold tongues. 

A widely accepted tool for studying entrainment is the 
phase response curve |2| . Specifically, the phase response 
function (infinitesimal phase response curve) indicates 
the phase shift of an oscillator due to an infinitesimal 
perturbation of a system variable [Jl . A classical problem 
in nonlinear dynamics uses phase response function and 
forcing waveform with which all important features of the 
entrainment process (e.g., locking range, defined as the 
width of the Arnold tongue at a given forcing amplitude) 
can be obtained for weakly perturbed systems [1]. 

Many applications require optimization of the entrain- 
ment process. This is often achieved by adjusting the 
forcing waveform to achieve a target entrainment fea- 
ture. A variety of control targets were explored: optimal 
input was determined for establishing fast entrainments 
[5|, circadian phase resetting p, |7|, starting/stopping of 
the oscillations [7|, |8|, and maximal resonancefenergy 
transfer) between the system and forcing signal (9| . Con- 
trol of deterministic p^ and stochastic |11L I12| neuronal 
spiking activity was achieved with phase modeling ap- 
proach combined with variational methods to optimize 
spiking time [10] and variance of firing rates |l2J . 

In this Letter, we give a full account for the inverse of 
the classical entrainment problem: what is the minimal 
power forcing waveform that produces efficient entrain- 



ment of a limit-cycle oscillator in weak forcing limita- 
tion? Although the quality of entrainment could involve 
features such as stability and basin of attraction, here we 
consider efficient entrainment as the occurrence of max- 
imum width (or minimum slopes) of the Arnold tongue. 
This 'locking range' quality marker has been commonly 
used for injection-lock micro integrated oscillators as well 
as phase- locked loop circuits [13]. We propose a versa- 
tile, efficient approach to obtain exact functional form of 
the optimal waveform provided that the response func- 
tion related to the forcing action had been established. 
The theoretically obtained optimal waveforms, which ex- 
hibit some unexpected symmetry relationships with the 
response function, are tested in a simple numerical model 
that include higher harmonics in the response function 
typically seen in strongly nonlinear oscillators. The ex- 
perimental applicability of the method is demonstrated 
with optimal entrainment of a chemical system, the os- 
cillatory electrodissolution of nickel in sulfuric acid. 

The entrainment process of a limit-c^le oscillator in 
weak forcing limit can be modeled by 
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^ = c^ + Z(V')/(m), 



(1) 



where ip is the phase of the oscillator, Z is the phase 
response function, and uo and Vt are the natural frequency 
of the oscillator and the forcing frequency, respectively. 
In weak forcing limit, Eq. ([1]) is further simplified by 
averaging [l5], as 



-^ = A. + rw, 



(2) 



where (p and Acj aregiven hy (j) = tp — Vtt^ and Aco' = 
rt — ix)^ respectively [l6|. The interaction function, r(0), 
is obtained from the forcing waveform / and the phase 
response function, Z, as T{(j)) = {Z {0 -\- (j)) f (0)) ^ where 
represents ftt and (•) denotes the average by 6 over 
its period 27r: (•) = (27r)~^^- d6. Entrainment occurs 
when the phase difference is locked, i.e., d(j)/dt = Ao; + 



r(0) = [TJ. The range of frequency difference, Acj, 
where solution for stable steady state exists for defines 
the locking range R[f] for a certain forcing waveform |17| . 
Therefore, the locking range is the difference between the 
maximum (at </) = ^+) and minimum (at == (j) ) values 
of r(^) where phase locked solution exists [18|. R[f] is 
thus given by r(0+) — T{(j)-). 

Now we are in position to formulate the optimal en- 
trainment problem mathematically: the optimal forcing 
waveform (/opt) maximizes the locking range R under 
certain constraints. A convenient practical constraint is 
the total power of the waveform over its period: {f{0)'^). 
Therefore, the optimal forcing waveform, /opt gives max- 
imal locking range for a given (constant) forcing power 
P. We consider this as a variational problem maximizing 
the functional form 



Here we test the above theoretical predictions by a 
simple model with the following response function. 



S[f]^R[f]-\{{f)-P), 



(3) 



where A is the Lagrange multiplier. Solution to the vari- 
ational problem, /h=, a suitable candidate for the optimal 
waveform, is obtained by ensuring that the first variation 
SS vanishes and the second variation S'^S is negative [l9|: 

M0) = {2X)-'{z{9 + 4>+) - z{e + 0_)}. (4) 

The Lagrange multiplier can be obtained by substituting 
the solution Eq. (J4]) in the constant power constraint 
{{fl) - P = 0): A = {l/2)^Q/P with Q ^ {{Z{e + 
0+) - Z{e + (t)-)Y)- Note that /* in Eq. g]) has zero 
average: (/*) = 0. 

To obtain /* of Eq. (|4]) the maximum (0+) and the 
minimum (</)_)(/) values of F with the forcing waveform 



r(0) = ./pjQ {z{e + 4>){z{e + 0+) - z{e + 0_)}> (5) 

have to be determined. The conditions for the maximum 
and minimum of F are as follows: 

F'((^±) = 0, F''(0+) < 0, and T"{(^_) > 0. (6) 

The first condition in Eqs. (|6]), combined with Eq. (|5]), 
gives 

{Z\0 + 0+)Z(^ + (/>_)) = {Z\e + A(/))Z(^)) = (7) 

where A^ = ^+ — ^_. (A0 is introduced to remove 
phase ambiguity). We shall refer to Eq. ([7j) as balanc- 
ing condition because this equation realizes optimality 
by balancing both terms in Eq. ^. The trivial A^ = 
solution to Eq. ([7j) is discarded because it does not allow 
entrainment (F((/)) = in Eq. ([5])). 

However, other solutions do exist in Eq. ([7j), because 

d{z'{e + A<p)z{e))/dA^\^^^, = -{z'{e)^) < 0, and 

{Z'{0 + A(j))Z{0)) is a periodic, bounded function of Acj) 
in a large class of systems [l|. In particular, we have 
found that in models with twice different iable, continu- 
ous Z a solution with A^ = n exists; we call this solution 
and the corresponding optimal waveform 'generic' |21|. 



Z{0) =sin(9 + asin(2(9). 



(8) 



This Z simulates the behavior of Stuart-Landau oscil- 
lator [1(1 with a = 0; therefore, we can consider a as a 
measure of distance from Hopf bifurcation that can in- 
troduce higher harmonics in the response function. The 
balancing condition of Eq. ([7]) for this model is explicitly 
written as [l + 4a^ cos(A^)] sin(A^) = 0. For \a\ < 1/2, 
there is only one (nontrivial) solution A(f) = it. The opti- 
mal waveform for this generic solution is independent of 
d'- fopt{0) = — \/2Psin^. Although the response function 
does contain second order harmonic for < \a\ < 1/2, 
this term does not affect the shape of the optimal wave- 
form. This finding suggests that with Z containing only 
weak second (or, in general even) harmonics the sinu- 
soidal forcing is the optimal since the generic solution 
to balancing condition always exists. For example, for 
systems close to Hopf bifurcation, which contain mostly 
first and weak second (even) harmonics in Z, the optimal 
waveform is sinusoidal. However, the odd harmonics do 
appear in the generic optimal waveform and thus systems 
with relatively strong third (and higher odd) harmonics 
in Z are expected to retain the odd harmonics in the opti- 
mal waveform. For the generic solution the locking range 
is calculated as i? = >/2P, which is again independent of 
a. 

In the range \a\ > 1/2, we have three different solu- 
tions for the balancing condition and thus three candi- 
dates for optimal waveform. The generic solution with 
A(j) = TT exists, however, there appear two additional, 
'non-generic' solutions satisfying cosA^ = — l/4a^. For 
these non-generic solutions the locking ranges are identi- 
cal: i? = (1 + 4a'^)VP/{2V2a) > VW. This shows that 
in the range \a\ > 1/2, the non-generic optimal wave- 
forms outperform the generic waveform; for large values 
of a, the improvement of locking range for the non-generic 
over generic waveforms increases approximately linearly. 
The non-generic waveforms are not purely sinusoidal and 
depend on the parameter a. Two non-generic waveforms 
for a = 0.95 are shown in Fig. 1(a). 

We have also verified these theoretical predictions, by 
using a standard genetic algorithm which numerically 
searches for /opt- Figure 1(a) shows all /opt obtained by 
the algorithm, for both cases of full (Eq. ([T])) and av- 
eraged (Eq. ([2])) phase models, with cj = 10, a = 0.95 
(> 1/2) and P = 0.5 or P = 10.0. The numerical al- 
gorithm found the exact same optimal waveform as pre- 
dicted by the theory, for both cases up to around P = 2.0. 
The numerical value of the locking range for these opti- 
mal solutions (maximum in R landscape in Figure 1(b)) is 
also the same as predicted by the theory; the non-generic 
optimal solution performs about 23.1% better than the 
simple sinusoidal (generic optimal) forcing, and also it 
performs about 89.8% better than the pulse forcing in 
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Figure 1: Optimal forcing /opt obtained in the model with the 
response function of Eq. (|8]). Panel (a) shows all /opt, after 
rescaling for comparison, obtained by the genetic algorithm 
for a — 0.95: green and blue curves for P — 0.5 with Eq. 
{[2]), dotted curves for P — 10.0 with Eq. ([!]). In this genetic 
algorithm, /opt is searched among all functions of the form 
cisin^ + C2sin(2^ + '02)- Panel (b) shows the landscape of R 
for a = 0.95 (> 1/2), with respect to (ci, '02). The peaks 1 
and 2 of the landscape correspond to the waveforms 1 and 2 
in panel (a), respectively. 



Fig. Wis)' Beyond P = 2.0, a small discrepancy appears 
in the optimal waveforms obtained with the genetic al- 
gorithm. However, the shape of bimodal landscapes of 
R[f] in Fig. [IJb) is preserved up to around P = 10.0 
for the case of Eq. ([1]), which suggests that, at least in 
this particular example, beyond the strictly weak forc- 
ing limit, the theoretical prediction of optimal waveform 
could be used as an initial candidate that can be further 
optimized with other techniques. 

The theoretical method for obtaining optimal forc- 
ing waveform was demonstrated in chemical experiments 
with Ni electrodissolution. The experiments were car- 
ried out in a standard three-electrode electrochemical 
apparatus with a 1 mm diameter Ni wire as working, 
a 1.57 mm diameter Ft coated Ti rod as counter, and 
Hg/Hg2S04/K2S04(sat) as reference electrode. The po- 
tential of the Ni wire was set to a circuit potential 
V = Vq -\- Afo{0), where Vq is a base potential, A is the 
amplitude of forcing, and fo{0) is the normalized forcing 
waveform with < /q >= 0.5. The current, proportional 
to dissolution rate was measured by the potentiostat. 
The frequency of the current oscillations was determined 
with the Hilbert transform method [1"| . The experiments 
were carried at 10 °C in 3 mol/L sulfuric acid solution. 
Further experimental details are given in previous stud- 
When a resistor of 1 kOhm is attached to the 
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Ni electrode, oscillations arise through Hopf bifurcation 
at Vq ^ 1.07 V due to the hidden negative differential 
resistance of the electrochemical system [24] . 

In accordance with previous experimental results [23| , 
at Vq = 1.100 V, slightly above the bifurcation point, 
the oscillations are smooth with a phase response curve 
that contains predominantly first harmonic components 
(see Figs. [2ja) and[2jb)). The experiments were carried 
out at fixed forcing frequencies and the amplitude of the 
waveform was increased until entrainment occurs where 
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Figure 2: Experiments: optimal entrainment of smooth os- 
cillators, Vq = 1.100 V. (a) Waveform, (b) Phase response 
curve (PRC). The PRC was obtained with pulsing the oscil- 
lations with an amplitude of -100 mV and width of 0.1 radian 
and measuring the phase advance ($) in 5 independent ex- 
periments. The error bars indicate standard deviations in the 
five experiments, (c) Three waveforms, sine-sawtooth-pulse, 
with which entrainments are tested, (d) Reduced entrainabil- 
ities of the three waveforms: P: pulse. Saw: sawtooth. Sin: 
sine. The detuning Acj/cj was set to ±5%; at each detuning 
two experiments were carried out. 



the critical amplitude, Ac, was recorded. To compare the 
forcing capability of different waveforms we determined 
the reduced entrainabilities as Er = \Auj\/{Acy^< /q >)• 
The average value of the reduced entrainability for fixed 
positive and negative detunings is proportional to the R 
value; therefore, it is a useful quantity to characterize 
the entraining capability of a waveform. In a typical set 
of experiments 5 % detuning was applied, i.e., Auj/uj = 
±0.05; it was confirmed that smaller detuning (2 %) gives 
Er values within 5 % relative error; experiments with 
larger detunings resulted in large changes of oscillation 
waveforms and thus the weak forcing assumption of the 
theory was not satisfied. 

We have tested three waveforms with smooth oscilla- 
tors; sine, sawtooth, and pulse waves are shown in Fig. 
[2](c) . The optimal waveform was determined by the above 
algorithm [25| ; only generic waveform exists and because 
the phase response curve has very small third and higher 
harmonics (less than 1%) we can consider the sine wave- 
form as optimal within experimental error. The entrain- 
abilities in Fig. [2](d) show that, as predicted, the sine 
waveform indeed outperforms the sawtooth and pulse 
waves. 

When Vq is increased to 1200 mV, the oscillations 
become moderately relaxational; the waveform and the 
phase response curve exhibits higher harmonics (see Figs. 
[31(a) and[3](b)). The analysis has revealed that three op- 
timal waveforms exist: one generic, and two non-generic 
waveforms given in Figs. [Sjc-e), respectively. In Fig. 
[Sje) we also show a waveform that is obtained by adding 
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that consider stability and basin of attraction, will be 
considered in a forthcoming publication. The proposed 
methodology provides a framework for efficient design of 
entrainment applications in electrical circuit technology 
(e.g., for injection- locked oscillators) as well as in biolog- 
ical pacemakers. 

I.Z.K. acknowledges financial support from Research 
Corporation Cottrel College Science Award. H.T. thanks 
Dr. Y. Ando and N. Sagayama for their help with genetic 
algorithms. 
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Figure 3: Experiments: optimal entrainment of moderately 
relaxational oscillators, Vo = 1.200 V. (a) Oscillation wave- 
form, (b) Phase response curve. The PRC was obtained simi- 
larly to that in Fig. [2l (c-e) Optimal waveforms, (c) Generic 
optimal waveform (Og). (d) Optimal waveform 1 (01). (e) 
Thick solid line: Optimal waveform 2 (02). Thin line: noisy 
optimal waveform 2 (02n). (f) Reduced entrainabilities of the 
tested waveforms. The ^rvalues are averages of experiments 
with ±2% and ±5% detunings. 



Gaussian random numbers with standard deviation of 0.1 
to the Fourier coefficients of the most optimal waveform; 
this waveform is referred to as noisy optimal 2 waveform. 
The reduced entrainabilities follow the theoretical pre- 
dictions; sawtooth, sine, and pulse waves are inferior to 
the entraining power of the optimal waveforms. As pre- 
dicted, the generic optimal waveform performs slightly 
worse than the non-generic optimal waveforms. The best 
waveform was optimal 2; note that the shape of this wave- 
form is strongly nontrivial and by the addition of small 
amount of noise to the waveform the E^ value decreases. 
Overall, optimal waveform obtained from the theory in- 
creased entrainabilities by about 50% to 90% relative to 
those obtained with sine/pulse/sawtooth waves. 

Construction of optimal waveform for entrainment was 
proposed and tested here with a single oscillator. The 
method, however, can be extended to a group of interact- 
ing oscillators where effects related to the collective phase 
response function [26] shall be considered. The optimal 
signal can also be applied in closed-loop feedback systems 
along with synchronization engineering [22| for seeking 
optimal target dynamics. A limitation of the method- 
ology is the requirement for weak forcing so that phase 
models can be applied. This limitation leads to an exten- 
sion of the method where the forcing signal is limited to 
small values. These extensions, along with other targets 
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